﻿#region License

// Copyright 2012-     Felix Obermaier
/*
 * This license governs use of the accompanying software. If you use the software, you
 * accept this license. If you do not accept the license, do not use the software.
 *
 * 1. Definitions
 * The terms "reproduce," "reproduction," "derivative works," and "distribution" have the
 * same meaning here as under U.S. copyright law.
 * A "contribution" is the original software, or any additions or changes to the software.
 * A "contributor" is any person that distributes its contribution under this license.
 * "Licensed patents" are a contributor's patent claims that read directly on its contribution.
 *
 * 2. Grant of Rights
 * (A) Copyright Grant- Subject to the terms of this license, including the license conditions
 *     and limitations in section 3, each contributor grants you a non-exclusive, worldwide,
 *     royalty-free copyright license to reproduce its contribution, prepare derivative works
 *     of its contribution, and distribute its contribution or any derivative works that you create.
 * (B) Patent Grant- Subject to the terms of this license, including the license conditions and
 *     limitations in section 3, each contributor grants you a non-exclusive, worldwide, royalty-free
 *     license under its licensed patents to make, have made, use, sell, offer for sale, import,
 *     and/or otherwise dispose of its contribution in the software or derivative works of the
 *     contribution in the software.
 *
 * 3. Conditions and Limitations
 * (A) No Trademark License- This license does not grant you rights to use any contributors'
 *     name, logo, or trademarks.
 * (B) If you bring a patent claim against any contributor over patents that you claim are
 *     infringed by the software, your patent license from such contributor to the software
 *     ends automatically.
 * (C) If you distribute any portion of the software, you must retain all copyright, patent,
 *     trademark, and attribution notices that are present in the software.
 * (D) If you distribute any portion of the software in source code form, you may do so only under
 *     this license by including a complete copy of this license with your distribution. If you
 *     distribute any portion of the software in compiled or object code form, you may only do so
 *     under a license that complies with this license.
 * (E) The software is licensed "as-is." You bear the risk of using it. The contributors give no
 *     express warranties, guarantees or conditions. You may have additional consumer rights under
 *     your local laws which this license cannot change. To the extent permitted under your local laws,
 *     the contributors exclude the implied warranties of merchantability, fitness for a particular
 *     purpose and non-infringement.
 *
 */

#endregion License

using System;
using NPack.Interfaces;

namespace NPack
{
    /// <summary>
    /// Cholesky Decomposition
    /// <para/>
    /// For a symmetric, positive definite matrix A, the Cholesky decomposition
    /// is an lower triangular matrix L so that A = L*L'.
    /// <para/>
    /// If the matrix is not symmetric or positive definite, the constructor
    /// returns a partial decomposition and sets internal flags that may
    /// be queried by the <see cref="PositiveDefinite"/> and <see cref="Symmetric"/> property.
    /// </summary>
    /// <typeparam name="TComponent">The type of the matrix components</typeparam>
    /// <typeparam name="TMatrix"></typeparam>
    public class CholeskyDecomposition<TComponent, TMatrix> : ICholeskyDecomposition<TComponent, TMatrix>
        where TComponent : struct, IEquatable<TComponent>, IComparable<TComponent>, IComputable<TComponent>, IConvertible, IFormattable
        where TMatrix : IMatrix<TComponent, TMatrix>, IEquatable<TMatrix>, IComparable<TMatrix>, IComputable<TMatrix>
    {
        private readonly IMatrixFactory<TComponent, TMatrix> _matrixFactory;
        private readonly TMatrix _inputMatrix;

        /* ------------------------
           Class variables
         * ------------------------ */

        /*
         * Array for internal storage of decomposition.
         */
        private readonly TMatrix _l;

        /*
         * Row and column dimension (square matrix).
         */
        private readonly int _n;

        /*
         * Symmetric and positive definite flag.
         */
        private readonly bool _isSym;
        private readonly bool _isPd;

        /* ------------------------
           Constructor
         * ------------------------ */

        /** Cholesky algorithm for symmetric and positive definite matrix.
        @param  A   Square, symmetric matrix.
        @return     Structure to access L and isspd flag.
        */

        /// <summary>
        /// Creates an instance of this class
        /// </summary>
        /// <param name="value"></param>
        /// <param name="matrixFactory"></param>
        public CholeskyDecomposition(TMatrix value, IMatrixFactory<TComponent, TMatrix> matrixFactory)
        {
            if (Equals(value, default(TMatrix)))
                throw new ArgumentNullException("value");

            if (matrixFactory == null)
                throw new ArgumentNullException("matrixFactory");
            _matrixFactory = matrixFactory;

            _inputMatrix = value.Clone();

            // Initialize.
            var a = value.Clone();
            _n = value.RowCount;
            _l = _matrixFactory.CreateMatrix(_n, _n);

            _isSym = (value.ColumnCount == _n);
            _isPd = true;

            var zero = value[0, 0].Zero;

            // Main loop.
            for (int j = 0; j < _n; j++)
            {
                var lrowj = _l.GetRow(j);
                var d = zero.Zero;
                for (var k = 0; k < j; k++)
                {
                    var lrowk = _l.GetRow(k);
                    var s = zero.Zero;
                    for (var i = 0; i < k; i++)
                    {
                        s = s.Add(lrowk[i].Multiply(lrowj[i]));
                    }
                    lrowj[k] = s = (a[j, k].Subtract(s)).Divide(_l[k, k]);
                    d = d.Add(s.Multiply(s));
                    _isPd &= (a[k, j].Equals(a[j, k]));
                }
                d = a[j, j].Subtract(d);
                _isPd &= d.GreaterThan(zero);
                _l[j, j] = Utility.Max(d, zero).Sqrt();
                for (var k = j + 1; k < _n; k++)
                {
                    _l[j, k] = zero.Zero;
                }
            }
        }

        /* ------------------------
           Temporary, experimental code.
         * ------------------------ *\

           \** Right Triangular Cholesky Decomposition.
           <P>
           For a symmetric, positive definite matrix A, the Right Cholesky
           decomposition is an upper triangular matrix R so that A = R'*R.
           This constructor computes R with the Fortran inspired column oriented
           algorithm used in LINPACK and MATLAB.  In Java, we suspect a row oriented,
           lower triangular decomposition is faster.  We have temporarily included
           this constructor here until timing experiments confirm this suspicion.
           *\

           \** Array for internal storage of right triangular decomposition. **\
           private transient double[][] R;

           \** Cholesky algorithm for symmetric and positive definite matrix.
           @param  A           Square, symmetric matrix.
           @param  rightflag   Actual value ignored.
           @return             Structure to access R and isspd flag.
           *\

           public CholeskyDecomposition (Matrix Arg, int rightflag) {
              // Initialize.
              double[][] A = Arg.getArray();
              n = Arg.getColumnDimension();
              R = new double[n][n];
              isspd = (Arg.getColumnDimension() == n);
              // Main loop.
              for (int j = 0; j < n; j++) {
                 double d = 0.0;
                 for (int k = 0; k < j; k++) {
                    double s = A[k][j];
                    for (int i = 0; i < k; i++) {
                       s = s - R[i][k]*R[i][j];
                    }
                    R[k][j] = s = s/R[k][k];
                    d = d + s*s;
                    isspd = isspd & (A[k][j] == A[j][k]);
                 }
                 d = A[j][j] - d;
                 isspd = isspd & (d > 0.0);
                 R[j][j] = Math.sqrt(Math.max(d,0.0));
                 for (int k = j+1; k < n; k++) {
                    R[k][j] = 0.0;
                 }
              }
           }

           \** Return upper triangular factor.
           @return     R
           *\

           public Matrix getR () {
              return new Matrix(R,n,n);
           }

        \* ------------------------
           End of temporary code.
         * ------------------------ */

        /* ------------------------
           Public Methods
         * ------------------------ */

        /// <summary>
        /// Gets whether the matrix symmetric and positive definite
        /// </summary>
        public Boolean PositiveDefinite
        {
            get { return _isPd; }
        }

        /// <summary>
        /// Solve A*X = B
        /// </summary>
        /// <param name="matrixB">A Matrix with as many rows as A and any number of columns.</param>
        /// <returns>X so that L*L'*X = B</returns>
        /// <exception cref="ArgumentException">Matrix row dimensions must agree.</exception>
        /// <exception cref="InvalidOperationException">Matrix is not symmetric positive definite.</exception>
        public TMatrix Solve(TMatrix matrixB)
        {
            if (matrixB.RowCount != _n)
            {
                throw new ArgumentException("Matrix row dimensions must agree.");
            }
            if (!_isPd || _isSym)
            {
                throw new InvalidOperationException("Matrix is not symmetric positive definite.");
            }

            // Copy right hand side.
            var nx = matrixB.ColumnCount;
            var x = matrixB.GetMatrix(0, _n, 0, nx);

            // Solve L*Y = B;
            for (var k = 0; k < _n; k++)
            {
                for (var j = 0; j < nx; j++)
                {
                    for (var i = 0; i < k; i++)
                    {
                        x[k, j] = x[k, j].Subtract(x[i, j].Multiply(_l[k, i]));
                    }
                    x[k, j] = x[k, j].Divide(_l[k, k]);
                }
            }

            // Solve L'*X = Y;
            for (var k = _n - 1; k >= 0; k--)
            {
                for (var j = 0; j < nx; j++)
                {
                    for (var i = k + 1; i < _n; i++)
                    {
                        x[k, j] = x[k, j].Subtract(x[i, j].Multiply(_l[i, k]));
                    }
                    x[k, j] = x[k, j].Divide(_l[k, k]);
                }
            }

            return x;
        }

        #region Implementation of IFactorization<TComponent,TMatrix>

        /// <summary>
        /// Gets the matrix which was factored (or decomposed).
        /// </summary>
        public TMatrix InputMatrix
        {
            get { return _inputMatrix.Clone(); }
        }

        #endregion Implementation of IFactorization<TComponent,TMatrix>

        #region Implementation of ICholeskyDecomposition<TComponent,TMatrix>

        /// <summary>
        /// Gets the left triangular factor
        /// </summary>
        public TMatrix LeftTriangularFactor
        {
            get { return _l.Clone(); }
        }

        /// <summary>
        /// Gets whether the L matrix is symmetric
        /// </summary>
        public bool Symmetric
        {
            get { return _isSym; }
        }

        #endregion Implementation of ICholeskyDecomposition<TComponent,TMatrix>
    }
}